clear
clear global variable
%载入网格数据
global mu own nei num d b S_p U_inlet c_s f_c V_c F_f P_out V
points=table2cell(readtable('points'));
faces=table2cell(readtable('faces'));
neighbour=readtable('neighbour');
owner=readtable('owner');
boundary=table2cell(readtable('boundary'));
p=points;
f=faces;
p(end-1:end,:)=[];
f(end-1:end,:)=[];
p=[str2num(char(strrep(p(:,1),'(',''))),cell2mat(p(:,2)),str2num(char(strrep(p(:,3),')','')))];
f=[str2num(char(strrep(f(:,1),'(',' '))),cell2mat(f(:,2)),cell2mat(f(:,3)),str2num(char(strrep(f(:,4),')','')))];
nei=cell2mat(table2cell(neighbour));
own=cell2mat(table2cell(owner));
nei(isnan(nei))=[];
own(isnan(own))=[];
own(2:end)=own(2:end)+1;
nei(2:end)=nei(2:end)+1;
f(:,2:end)=f(:,2:end)+1;
[f_c,c_s]=face_centered(f,p);%求面心及面向量
[Array,c_s_]=arrayBuild(c_s);%建立稀疏矩阵，使用稀疏矩阵表示面向量
[V_c,V,Cell]=volume;%求体心与体积及单元对应的面

%设置初始值及计算时间
U_inlet=[1,0,0];%输入初始值
P_out=0;
rho=1;%输入密度
mu=0.01/rho;%输入动力粘度
P.t0=zeros(num,1)/rho;%各网格初始时刻压力值
U.t0=kron(ones(num,1),U_inlet);%各网格初始时刻速度场
delta_t=10^-4;
t_end=10^-4;

%建立边界条件
num_=str2double(cell2mat(boundary(1)));
bt=find(strcmp(boundary, '{' ));
nfaces=find(contains(boundary ,'nFaces'));
for i=1:num_
    field=cell2mat(boundary(bt(i)-1));
    b.(field)=[str2double(cell2mat(strrep(boundary(nfaces(i)),'nFaces',''))),str2double(cell2mat(strrep(boundary(nfaces(i)+1),'startFace','')))+1];
end

%进入迭代
n=t_end/delta_t;
d_err=zeros(n,1);
for i=1:n
    t=delta_t*i;
    dt=sprintf('t%d',i);%预测步
    dt_=sprintf('t%d',i-1);%前一步
    U.(dt)=zeros(num,3);
    %构建该时刻各参数
    [F_f,F_f_c]=face_flux(U.(dt_));%面通量,中心线性格式
    S_p=U.(dt_)/delta_t;
    A=AarrayBuild(F_f_c,V,delta_t,c_s_);    
    %预测速度
    p_A=PBuild(P.(dt_));
    U.(dt)=A\(S_p-p_A);
    d_=ones(num,3);
    d_err(i)=max(d_);
    error=10^-5;
    j=0;
    while j<=3
        j=j+1;
        Ap=diag(A);
        An=A-diag(Ap);
        HbyAp=(-An*U.(dt)+S_p)./(Ap*ones(1,3));
        P.(dt)=PLaplacesolver(U.(dt),HbyAp,Ap,c_s_);
        p_A=PBuild(P.(dt));
        p_A=p_A./(Ap*ones(1,3));
        Ud=HbyAp-p_A;
        d_=abs(Ud-U.(dt));
        U.(dt)=Ud;
        d_err(i)=max(d_);
        if d_err(i)<=error
            break;
        end
    end
end       
